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Myelin figures are densely packed stacks of coaxial cylindrical bilayers that are unstable to the for- 
mation of coils or double helices. These myelin figures appear to have no intrinsic chirality. We show 
that such cylindrical membrane stacks can develop an instability when they acquire a spontaneous 
curvature or when the equilibrium distance between membranes is decreased. This instability breaks 
the chiral symmetry of the stack and may result in coiling. A unilamellar cylindrical vesicle, on the 
other hand, will develop an axisymmetric instability, possibly related to the pearling instability. 

PACS numbers: 87.16.Dg 



I. INTRODUCTION 

Amphiphilic molecules in water self-assemble into a va- 
riety of structures, including micelles, bilayer vesicles, 
and stacks of bilayer membranes. When a dehydrated 
lump of amphiphile is brought into contact with water, it 
will develop a large number of tubular structures which 
grow to lengths many times larger than their radii [Q. 
These tubes, called myelin figures, are densely-packed, 
nested stacks of coaxial cylindrical bilayers [|| . Typically 
the outer radius of such a stack is on the order of tens 
of microns and the inner radius can be as small as or 
smaller than 0.2/im Q (see FIG. [I]). Myelin figures tend 
to be straight over lengths many times longer than their 
radii, but eventually bend or fold back on themselves || . 
Though they appear to have no intrinsic chirality, they 
are observed to form coils or double helices ||] • 

Coiling has also been observed in cylindrical stacks of 
binary mixtures j(| , and egg- yolk phosphatidylcholine . 
In the first system, Ca 2+ ions were added to myelin fig- 
ures formed from a binary mixture of cardiolipin and 
phosphatidylcholine. Above a critical concentration of 
ions, tightly-packed, single helices were formed. Coiling 
was also observed when the system was hydrated with a 
solution of Ca 2+ ions. In the latter system, the myelin 
figures would be straight initially but form helices and 
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FIG. 1: A myelin figure is a multilamellar cylindrical stack 
of bilayers. Typical radii are on the order of 10 /im. 
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other structures after they had been growing for several 
minutes. It was presumed that the observed coiling in 
both systems was due to surface adhesion. However, coil- 
ing has also been observed recently in another system of 
myelin figures hydrated with a polymer solution in which 
it was demonstrated that surface adhesion played a neg- 
ligiable role in the coiling mechanism [|| ||. Coiling in 
this system was attributed to an interaction between the 
polymer and the membranes which resulted in a local 
spontaneous curvature. In a model was developed 
in which the polymer is allowed to diffuse freely along 
the bilayer surfaces. Once the polymer concentration is 
above a critical concentration, the straight myelin figure 
becomes unstable to one which is maximally coiled. 

Unilamellar tubes in the presence of the same polymer, 
on the other hand, are known to exhibit a pearling insta- 
bility 0. Here, the tubes are formed first and a drop of 
polymer solution is added later. In the early stages of the 
instability, the tubes develop a periodic oscillation along 
their length which leads, eventually, to the formation of a 
string of spheres. A different pearling instability has also 
been observed in unilamellar tubes excited with optical 
tweezers J || and has been attributed to induced surface 
tension |9, and unilamellar tubes with oil incorpo- 
rated into the bilayer [ pd| . 

In this paper, we calculate the bending and interac- 
tion energies for undulations of a cylindrical multilamel- 
lar stack. We then study the instabilities generated by 
perturbing the stack in various ways: (1) we add sponta- 
neous curvature to all the bilayers, and (2) we decrease 
the interaction equilibrium distance between the layers 
after the stacks have been formed. We show that the 
most unstable mode (neglecting dynamics) breaks the 
chiral symmetry of the cylinders and introduces a twist- 
ing undulation, which may be related to coiling. In con- 
trast, spontaneous curvature on a unilamellar cylindrical 
vesicle results in periodic oscillation along the cylinder 
axis, as in the pearling instability. 

A similar instability in a smectic-A liquid crystal 
confined within a cylindrical channel with homeotropic 
boundary conditions was studied in Q. The bound- 
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ary conditions cause the smectic layers to arrange them- 
selves into co- axial cylinders much like myelin figures. 
The presence of the confining walls changes the nature of 
the instability, however, since the outer layer of a myelin 
figure is not fixed to be a cylinder. Also, those authors 
consider only undulations which vary along the length of 
the cylinder which explicitely prohibits instabilities that 
break chiral symmetry. 

The remainder of the paper is organized as follows. In 
section II we calculate the total energy of a cylindrical 
multilamellar stack to second order in arbitrary undula- 
tions of the radius. In section III we demonstrate the 
existence of the instabilities. In section IV, we discuss 
the nature of the instabilities. Finally, we summarize in 
section V. 



II. ENERGY OF A CYLINDRICAL STACK 

A. Geometry 

We will assume the bilayers in the cylindrical stack to 
be equally spaced jL9|, though the actual distribution of 
layers is unknown. At length scales larger than molecular 
lengths, it is appropriate to model a bilayer membrane 
as a two dimensional surface. There are two equivalent 
representations of the membrane surfaces that we will 
use. 

On one hand, we can introduce local coordinates on 
the cylinders (9, z) and a vector function X(9, z) which 
gives the position in space of any point on the membrane. 
The n th membrane in the stack can then be described by 
X(9,z) = r n (9, z)(cos(9)x+sin(9)y)+zz where r n (6, z) = 
ro is independent of 9 and z for a cylinder of radius tq . 

We generalize this to r n (9,z) = an + h n (9,z) where 
a is the spacing between the layers, and h n (9,z) is an 
undulation of the cylinder radius. Since the membranes 
in the stack are very close together, the undulations of 
the layers are severely constrained and h n (9,z) should 
be slowly varying radially. Therefore, to simplify our 
calculations, we specialize to the case where h n (9,z) is 
independent of n. 

The induced metric, g a p = 9,z), is defined 

by g a p = d a X ■ dpi where d e = ■§§ and d z = ^. 
The differential surface area element is dA = dBdz^fg 
where g = det g a p- The mean curvature and Gaus- 
sian curvature are defined to be the trace and determi- 
nant respectively of the curvature tensor, K a p, given by, 
K a ph = D a dpX where D a is the covariant derivative. 
Then H = |AA| where A = g a ^D a dp is the Laplace- 
Bcltrami operator acting on a scalar. In two dimensions, 
J cL4det K a p is a topological invariant. 

An alternative way to represent the layers is with a 
function (f>(r, 9, z) — r — u(r, 9, z). The n th membrane is 
then given by the solution to the equation cj>(r, 9, z) = a$n 
for some ao. It is straightforward to show that u(r, 9, z) = 
(a-a o )r+a o h(0,z) ^ggQj-j^gg a cylinder with a radius given 



by r n (9,z) = an + h{9,z). In this representation, we 
write the unit normal as n ~ j^y- This gives the unit 
normal as a function of (r, 9, z) so we use h(r n (9, z), 9, z) 
for the unit normal on the n th membrane. Finally, we 
write H(r,9,z) = V.n and H n {9,z) = H(r n (9 : z),0, z) 
gives the mean curvature on the n th membrane. This 
representation is similar to that used for smectic liquid 
crystals (see, for example, Q). 

B. Bending Energy 

The theory of the bending energy of amphiphilic mem- 
branes is well-established for small curvatures |l4|]. The 
bending energy is given by, 

Eb = \\ dA{ - H ~ H °^ 2 + R J dAK ^ 

where H is the mean curvature, K — det K a p is the 
Gaussian curvature, and Hq is the spontaneous curva- 
ture. Spontaneous curvature can arise from an asymme- 
try between the two monolayers of the membrane. Here 
we will assume that H = 0. Since J dAK is a topolog- 
ical invariant and we do not allow topology change, we 
suppress this term. 

We assume that the surface area per molecule is con- 
stant. Though myelin figures are nonequilibrium struc- 
tures, we will also assume that their rate of growth, 
— 0.3/ito/s ||, is slow enough to allow us to treat them 
as being in equilibrium at any given time. These assump- 
tions together give us the constraint that the total surface 
area of the layers should remain constant. To impliment 
this constraint, we add an effective surface tension of the 
form a J dA where J dA is the surface area of the mem- 
brane and the surface tension, er, is actually a Lagrange 
multiplier chosen to enforce the constraint. 

This single constraint, however, is not enough to en- 
sure that a cylindrical membrane is an extremum of the 
bending energy. If we neglect the endcap energy (which 
is reasonable for long cylinders), the mean curvature for 
a cylinder of radius p is 1/ p. Thus, the bending energy 
of a cylinder can always be lowered by increasing the ra- 
dius, and we can always do so while conserving area by 
decreasing the length. An additional contribution to the 
free energy is necessary for a cylindrical membrane to be 
stable. 

Since the mechanism which stabilizes a myelin figure 
is still poorly understood, there are a number of possi- 
ble approaches we can take to this problem. On the one 
hand, Sakurai and Kawamura jj] suggest that a cylin- 
drical stack can be stabilized by the interaction energy, 
which forbids the layers from increasing their separation. 
This still allows for the possibility of increasing all the 
radii without changing their separations. A second pos- 
sibility is that of osmotic pressure which keeps the volume 
between the layers constant. In |ll|, the authors use os- 
motic pressure to stabilize a cylindrical vesicle and their 
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results have a simple generalization to cylindrical stacks. 
However, for mathematical simplicity, we choose to solve 
this by adding a linear tension (of the form —tL). By 
choosing the linear tension and surface tension properly, 
we can ensure that a cylindrical stack of any length and 
distribution of layers is an extremum of the energy. In 
section V, we comment on how conserving the volume 
between the layers changes our results. 

Finally, the effective energy of each membrane layer, 
without interactions between the layers and with Hq = 0, 
is given by, 



e B = a J d 2 Xy/g + ^ J d 2 Xy/gH 2 - tL 



(2) 



It will be convenient to use the Fourier series represen- 
tation for the function h{9, z). Since the cylinder length 
is much longer than the radius, we will neglect the con- 
tribution of the ends of the stack to the energy and will 
impose periodic boundary conditions along the cylinder 
axis. This allows us to write h(9, z) = J2 m j b m je me+lqjZ 
where qj = -jf-, L is the cylinder length, and both m 
and j are integers. The reality of h(9, z) imposes the 
condition b* m j = b- m -j on the Fourier coefficients where 
* signifies the complex conjugate. The coefficient &oo rep- 
resents changes of radius of all the bilayers. This coef- 
ficient is redundant with the layer separation a and the 
core radius. 

We calculate a and t for each layer by imposing the 
conditions that the first derivative of the bending energy 
of a straight cylinder with respect to both radius and 
length should vanish. This gives us a — and t = 



where p is the radius of the cylinder. Expanding 
the bending energy to second order in h{9, z) and its 
derivatives and using these values for a and t results in 
the expression, 
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for the bending energy of each bilayer, where p = an is 
the n th membrane radius. 

Notice that the first order variation of the bending en- 
ergy with respect to the undulations h{9, z) vanishes since 
there are no terms in the bending energy linear in b m j. 
Therefore, we can conclude that a straight cylinder is an 
extremum of this bending energy. Also, notice that the 
undulation mode with j = and m = 1 has zero bend- 
ing energy. Since the shape change associated with this 
mode causes a change in the mean curvature of higher 
than second order, it is typically associated with trans- 
lation of the cylinder normal to the long axis. 

Finally, to find the total bending energy we sum over 
the membranes. Since they are very dense, we go over to 
the continuum limit and replace with - J R dp where 



R is the radius of the stack, and r c is the radius of the 
core. The final result for the bending energy is, 



E B = 

R 



,j ( TO 2 - l) 2 

2r? 



2 q] + 2m 2 q 2 ln(R/r c 
where we have used the fact that r c <C R. 
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C. Interaction Energy 

The largest contribution to the interaction energy 
comes from nearest-neighbor membranes. Consider the 
shortest curve between the two membranes which is nor- 
mal to both surfaces. If the membrane separation and 
size of undulations are small, this curve is nearly straight 
and, for any point (9, z) on one membrane with normal 
curve passing through it, the interaction energy gets its 
dominant contribution from the point on the other mem- 
brane which lies along the curve. Let a be the length of 
the normal curve between the two membranes, a be the 
separation of the layers with no undulations, and a be 
the equilibrium separation of the layers. For simplicity, 
we will assume that a and a are constant (in particular, 
they are independent of n) . 

To find d{9, z), we represent the membrane surfaces as 
solutions to the equation a^n = 0(r, 9,z) = r — u(r, 9, z). 
Then, noting that § c in = in the absence of defects 
and that n = we can write a(#, z) — k, a[l — 

(f-K^) 2 )] 
We can then expand the interaction energy in powers of 

(a — o,q) /ao to second order. The first order term vanishes 
since ao is the equilibrium separation, so the interaction 
energy density is -f [§7 — ^(Vu) 2 ] 2 where B is the bulk 
modulus. This must be integrated over the surface of 
each membrane in the layer. This interaction energy is 
similar to the familiar de Gennes elastic term for the 
free energy of a smectic liquid crystal p"3| . However, our 
bulk modulus B has units of energy per area rather than 
energy per volume as it does in a liquid crystal. 

To properly choose a and t, we should take into ac- 
count the interaction energy of the membranes. Do- 
ing so gives a correction to our value for a at nonzero 
7 = 0=2* . In fact, this gives a = - f ^g+ffl' 

where p = (1 + j)aon is the radius of the n th membrane. 

Using our value for u(r, 9, z), the resulting interaction 
energy per layer is, 
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7/2)$> m ,f[- 



(5) 



If there are no undulations, the interaction energy is 
zero even for nonzero 7. This comes about because of the 
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correction to a. Without this correction, a cylindrical 
stack with a ^ a will not be an extremum of the energy. 

We finally note that if 7 = 0, there is no contribution 
from the interaction energy to second order. This is be- 
cause we are not allowing the undulation of the radius, 
h(9, z), to vary from bilayer to bilayer so, to second order, 
they remain equidistant. 



III. INSTABILITIES 



A. Spontaneous Curvature 



Adding the terms, 



E c = -k / d'x^fgHHo 



d 2 x^H 2 (6) 



to the bending energy of each layer gives each membrane 
a spontaneous curvature, Hq. In order to conserve sur- 
face area the surface tension, a, will adjust itself to can- 
cel the second term exactly. Without changing the linear 
tension, however, there will appear terms in the total 
energy that are linear in the b m j arising from the first 
term of equation (|^). This corresponds to an instability 
that decreases the radii of the lamellae. However, the in- 
teraction energy of the stack requires that the change in 
radius of the layers be independent of the radii of the lay- 
ers. Conservation of surface area implies that the change 
in length, SL, be related to the change in radius, SR, by 
the expression SL = —SRL/p. 

Since 5R < 0, we see that the length of each layer is in- 
creased. However, the change of length is smaller for the 
layers with larger radii. The endcap of the inner layers 
will, therefore, be able to reach the endcaps of the outer- 
layers very early in this instability. This will prevent any 
large changes of length in the cylindrical stack. Fixing 
the length of the cylindrical stack can be accomplished 
by adjusting the linear tension to be t = 2j ^ L — 2-kkHq. 
With these changes, terms that are first order in the b m j 
vanish. Expanding the resulting energy to second order 
in the undulations results in the following additional term 
to equation (||), 

Cm 2 



e c = —2ttL 



I bmj |' 



with C = and D = kHq. Approximating the sum over 
membranes with an integral, the total stack energy is 
calculated to be ^fr^ = Ylmj \b m j\ 2 E m j, where, 



Ernj — 



-Cm 2 



K(m 2 



Dq 2 R 
-I) 2 



4r-2 



nq^R 2 



Km <7 7 ln(i?/r c 



(8) 



Since L is large, qj is well-approximated as a continu- 
ous variable, q. For m = 0, we have, 



0q 



K 
4f2 



(9) 



Minimizing this to find the most unstable mode gives 
2 _ 2H^ with an onse t threshold of H > ^— . For m = 1, 



1 R 

we have, 



E lq = —kHoqR 



Kq i R 2 



Kq 2 \n(R/r 



(10) 



This expression is negative whenever q 2 < q 2 



4H 

r 



■^\n(R/r c ) and H > ^ln(R/r c ). This threshold is 
much smaller than the threshold for the m = insta- 
bility. Finally, larger values of m can be analyzed like 
the m = case, yielding even larger thresholds for in- 
stability. Notice that q 2 is proportional to H — H 0c , so 
the critical wavelength is very long for spontaneous cur- 
vatures just above H$ c . The most unstable wavelength 
can be found by minimizing E\ q . Doing so gives the most 
unstable mode to be q 2 = 2^ - 2R 2 \n{R/r c ). 

We can do a similar analysis on a unilamellar mem- 
brane by using the energy for a single layer. The pro- 
cedure above establishes the existence of an instability 
when Hq > 1/p where p is the cylinder radius for both 
m = and m = 1. The m = 1 instability starts with 
long wavelength modes, however, while the m = insta- 
bility has an onset wavelength of 2irp and is expected to 
dominate. 

This result qualitatively agrees with Ou-Yang and 
Helfrich's calculation of instabilities in cylindrical vesi- 
cles M] stabilized by osmotic pressure. One crucial dif- 
ference, however, is that the m = 1 mode never goes 
unstable in their calculation. The discrepancy must be 
in our choice of linear tension over osmotic pressure to 
stabilize the cylindrical vesicle. 



B. Decreasing the interaction equilibrium distance 

We now consider the situation when the equilibrium 
distance between the membranes is decreased. This re- 
quires special care because the mode which relaxes the 
equilibrium distance by reducing the radii uniformly and 
increasing the length has been neglected by assuming 
h(9, z) is independent of radius. This relaxation mode 
carries a linear dependence on the radius, p, and will 
need to be considered explicitely. 

There are two effects which compete against the com- 
plete relaxation of the interaction energy through this 
mode. The first is the curvature energy which resists de- 
creasing the radii of any of the membrane layers. The sec- 
ond is the endcap energy which resists changes of length. 
The stack will be unable to completely relax its inter- 
action energy because of these effects, and the resulting 
separation between the layers, a, will be larger than the 
new equilibrium distance, do- This is mathematically 
equivalent to a swelled stack, which for a flat, lamellar 
phase results in a Helfrich-Hurault type instability |lq| . 

For a cylindrical stack, this results in an additional 
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contribution to equation (g) of, 
ej = -2ttL^|6. 



P 



(11) 



where C = D = f ^g±^ 



Summing over all 
membranes in the stack gives a total energy i? = 



^E^W^, where, 



Dm 2 ln(i?/r c ) + — q 2 i? 2 -r - 
2 

+ —r<& + K"i 2 9? ln(-R/r c ) 



«(m 2 — l) 2 



4r 2 



(121 



We again study the stability of the energy for different 
choices of m while taking qj to be continuous. For m = 0, 
we have, 



Eoq — 



Dq 2 R 2 



Ar 2 



nR 2 q 



2„4 



(13) 



This equation gives a most unstable mode of q 2 



D 



and an instability whenever D > 
have, 



Br, 



For m = 1, 



Si, - -Dln(R/r c 



Dq 2 R 2 kR q 



2„4 



which is negative when q 2 < q 2 



nq 2 ln(R/r c ) 
(14) 



2 4 
Thus, the reality of q requires that 



d 2 , iMflWF 
7?" H 

£> > £> c = 0, or 7 > 0. Again, this threshold is smaller 
than the m = threshold and the thresholds for in- 
stability for m > 1 are even higher than the m = 
case. We see then that, as D is increased, the first mode 
to destabilize is the one with m — 1 and q 2 ss for 
small 7. We find the most unstable wavenumber to be 
q =R-^\n(R/r c ). 

The critical wavelength, 2ir/q c — (27tAv2)/a/7 where 
A = y/B/n is the smectic penetration length, can be 
quite short. A 0.1 percent reduction in equilibrium dis- 
tance results in a critical wavelength of approximately 
280 A. If A is on the order of the layer spacing and there 
are several hundred layers in the cylindrical stack, the 
critical wavelength will be on the order of a tube radius. 
A one percent reduction in equilibrium distance results 
in a critical wavelength smaller than a single tube radius. 



IV. COILING OF MYELIN FIGURES 

Whenever q > 0, there are two types of behavior that 
can result, which we refer to as twisting (or coiling) and 
writhing. Twisting can be represented by a model un- 
dulation of the form h(9,z) — acos(6 + qz), for ex- 
ample, which breaks the chiral symmetry of the cylin- 
ders. Writhing, on the other hand, can be represented 





(a) 



(b) 



FIG. 2: (a) A three dimensional representation of a twist- 
ing undulation, (b) A three dimensional representation of a 
writhing undulation. 



by h{9, z) — acos(#) cos{qz) (see FIG. ||). Fourier coeffi- 
cients which give twisting for m = 1 and some choice of 
q are b\ q = a/2 and = with the remaining coeffi- 
cients zero or set by the reality condition on h(9,z). In 
the case of writhing, we have bi q = a/A and = a/4. 
The normalization factors of 2 and 4 are necessary so 
that the amplitudes of both undulations are equal to a. 
Now we can write the energy for both instabilities in the 
form, 



E,, 





a 2 




2 






— Emq 


a: 2 


4 



(15) 
(16) 



So we see that whenever E mq < 0, the most unstable 
case is twisting and not writhing. Since the m = 1 mode 
moves the center axis of the cylinders in the stack, we 
speculate that this mode may result in coiling as seen in 
experiments. 

In both types of instabilities discussed, the first mode 
to go unstable is the mode with m = 1 and q — 0. We 
expect that, when dynamics is taken into account, this 
q = mode will grow more slowly than all other wave- 
lengths. Thus, the dominant wavelength will be at some 
q > 0. It is not clear, however, how dynamics will af- 
fect whether the stack twists or writhes, and whether the 
behavior will lead to coiling. 
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V. SUMMARY AND DISCUSSION 

In this paper, we considered a nested stack of cylindri- 
cal, interacting membranes. We studied the instabilities 
generated by perturbing the stack in two ways: (1) we 
add spontaneous curvature to all the bilayers, and (2) 
we decrease the interaction equilibrium distance between 
the layers after the stacks have been formed. In both 
cases, we find that the first modes to go unstable are 
long wavelength modes which break the rotational sym- 
metry of the cylinder. We find that the most unstable 
modes break chiral symmetry and may be related to the 
coiling which is observed in myelin figures. 

We should stress that we have studied only two in a 
larger class of instabilities. If we imagine perturbing the 
system by adding a term to the free energy that prefers 
local undulations of the layers, we are constrained, by 
dimensional analysis, to write down, 

E P = 2ttL^ |6 mj f(-C7mV - Dq)p^ 2 ) (17) 
»'./ 

to quadratic order in b m j, m, and qj. The units of C and 
D (assumed the same) determine the value of v. As we 
have shown, a chiral symmetry breaking instability will 
occur when v = — 1 or y = - 2. Different mechanisms 
might then give rise to similar instabilities, depending 
only on the dimensions of their coefficients and the gen- 
eral condition that they result in local undulations. We 
expect, therefore, that there may be many instabilities of 
nested cylindrical stacks that lead to similar behavior. 

The instability resulting from decreasing the layer 
equilibrium distance bears some resemblence to an in- 
stability due to deswelling a myelin figure discovered by 
Chen, et al. jl7| |l8). They consider myelin figures whose 
outer layers only have been deswelled. In this case, the 
in-plane elasticity of the layers becomes relevant. They 



find that deswelling a myelin figure in this way results in 
an axisymmetric instability (m = 0). 

It is important to note that allowing 7 or Hq to vary 
with radius can change the results completely. We see 
this already if we consider adding spontaneous curvature 
only to the outermost layers of the stack. In that case, 
the additional energy will be approximately proportional 
to the energy of a cylindrical vesicle. As we have already 
seen, an energy of this form has a short wavelength ax- 
isymmetric instability (m = 0) and a long wavelength 
helical instability (m = 1), both at the same threshold. 

We have already mentioned that our results for a 
spontaneous curvature-induced instability in a cylindri- 
cal vesicle disagree with those of jll| , which uses osmotic 
pressure rather than linear tension to stabilize the vesi- 
cles. This sensitivity to our choice of constraints carries 
over to a cylindrical stack, also. A cylindrical stack with 
fixed volume between the layers which acquires sponta- 
neous curvature on all of its layers will not develop an 
instability at m = 1 as in the case of a cylindrical vesicle. 
The first mode to destabilize must then be m = 0. We 
have checked, however, that using osmotic pressure on 
a cylindrical stack whose equilibrium distance has been 
decreased still results in an instability with m — 1 and 
onset 7 > 0, but does change the value of q c . 
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